Inferring direction of associations between histone modifications using a neural processes-based framework

Summary Current technologies do not allow predicting interactions between histone post-translational modifications (HPTMs) at a system-level. We describe a computational framework, imputation-followed-by-inference, to predict directed association between two HPTMs using EpiTOF, a mass cytometry-based platform that allows profiling multiple HPTMs at a single-cell resolution. Using EpiTOF profiles of >55,000,000 peripheral mononuclear blood cells from 158 healthy human subjects, we show that neural processes (NP) have significantly higher accuracy than linear regression and k-nearest neighbors models to impute the abundance of an HPTM. Next, we infer the direction of association to show we recapitulate known HPTM associations and identify several previously unidentified ones in healthy individuals. Using this framework in an influenza vaccine cohort, we identify changes in associations between 6 pairs of HPTMs 30 days following vaccination, of which several have been shown to be involved in innate memory. These results demonstrate the utility of our framework in identifying directed interactions between HPTMs.


INTRODUCTION
Histone post-translational modifications (HPTMs) play a vital role in the regulation of gene expression, cell differentiation, and different processes centered around DNA. 1 Dysregulation of HPTMs has been implicated in human diseases such as cancer, 2,3 infectious diseases, 4 mental illnesses, 5 and autoimmune disorders. 6,7 HPTMs are also known to play an important role in the immune response following vaccination. 8,9 Interaction of multiple HPTMs has been recognized as the language of histone crosstalk. 10,11 Despite the advances in understanding the chromatin organization and transcriptional regulation through HPTMs, histone crosstalk has been understudied due to technological limitations of ChIP-seq 12 to investigate more than 2-3 HPTMs at once. iScience Article approaches have been used to impute multi-panel cytometry data. For instance, CyTOFmerge uses a k-nearest neighbors (kNN)-based approach to create an imputation function for each subject in the dataset. 16 Such subject-specific models are advantageous since individuals are exposed to different environments that are known to have different effects on epigenetics and HPTMs. 13,17 However, it is difficult to interpret the importance of input variables in kNN because it is non-parametric. Other approaches use parametric non-linear models such as support vector machines, boosting trees, and canonical neural networks. 18 Although these parametric models can be interpreted to quantify the importance of each input variable, they use the same imputation function for all subjects and do not distinguish between them, which may not capture the effects of difference in environmental exposure between individuals.
Neural process (NP), a recently described neural network-based model, combines the advantages of both kNN and parametric non-linear models. 19,20 Similar to kNN-based models, NP-based models use subjectspecific functions for imputation. Because NP-based models are parametric, the importance of each input variable can be interpreted. Importantly, NP models are scalable to complex functions and large datasets. 19 Because the biological processes regulating the interactions between HPTMs are complex, we hypothesized that NPs would more accurately impute HPTMs and CPMs, and reveal biologically meaningful relationships between them.
Here, we adapted the neural network architecture from NP to impute multi-panel mass cytometry-based HPTM abundances using subject-specific parametric models. Using more than 55 million human peripheral mononuclear blood cells (PBMCs) from 158 healthy subjects across 25 independent experiments, we evaluated linear regression (LR)-, kNN-and NP-based methods for accuracy of imputation. We also developed an interpretation framework to infer the direction of association between two HPTMs through systematic perturbations. We reproduced several known HPTM associations and identified several previously unidentified associations between pairs of HPTMs. Most of the associations were conserved across all healthy individuals. By combining the inferred pair-wise associations, we created a system-wide directed network of HPTM associations from the NP models. Finally, we demonstrated the utility of the NP models and the interpretation framework by identifying HPTMs whose associations were modified in response to the trivalent inactivated seasonal influenza vaccine (TIV) in healthy adults.

Data and model types
We profiled 55.6 million PBMCs from 158 healthy subjects (13-80 years old) across 25 EpiTOF experiments and measured 38 HPTMs and histone variants, 2 core histones (H3 and H4), and 11 cell phenotypic markers (CPMs) across two panels referred to as the methylation panel and the acetylation panel (Figures 1A and S1;  Table S1, STAR Methods). We have previously described extensive validation of antibodies for HPTMs using Western blot, flow cytometry, and different cell lines with in vitro manipulation. 13 CPMs are expressed on the cell surface and are used to determine the immune cell sub-type, whereas HPTMs are inside the cell nucleus. Given these biological differences between CPMs and HPTMs, we compared the accuracy of LR, kNN, and NP models for imputing a CPM or an HPTM using three imputation tasks: impute an HPTM using CPMs (task 1), impute a CPM using HPTMs (task 2), and impute an HPTM using other HPTMs on the same panel (task 3) ( Figure 1B). We used the training and validation sets for creating imputation models. Next, we locked each model and evaluated its accuracy using the test set. We evaluated 98 combinations of inputs and outputs for 294 models across the three tasks for each of the ML methods (LR, kNN, and NP; Figure 1C and Table S3). We assessed the accuracy of a model using the coefficient of determination (R 2 ), where an R 2 of 1 indicates an accurate imputation whereas an R 2 of 0 indicates a failure to impute accurately. Importantly, we only report results from the test set experiments because the imputation models are expected to have high accuracy in training and validation set experiments.
Cell phenotypic mark and histone post-translational modification abundances are mostly independent of each other When imputing an HPTM using CPMs in task 1, each of the three ML methods had low to moderate accuracy. Specifically, the mean R 2 for LR, kNN, and NP models were 0.19 (range: 0-0.4), 0.28 (range: 0-0.55), and 0.29 (range: 0.01-0.56), respectively ( Figure 1C). Although NP had significantly higher mean R 2 compared to LR and kNN models (p < 0.002), the low to moderate imputation accuracy suggests that CPMs alone are not sufficient to impute HPTMs across panels accurately.
Similarly, when imputing a CPM using HPTMs in task 2, each of the three methods had low to moderate accuracy. Specifically, the mean R 2 for LR, kNN, and NP were 0.15 (range: 0.02-0.47), 0.28 (range: 0.10-0.57), and 0.30 (range: 0.08-0.60), respectively ( Figure 1D). NP had a significantly higher mean of R 2 compared to LR and kNN models (p < 0.09). Taken together, our results indicate that the abundances of CPMs and HPTMs are largely independent of each other. These results are in line with single-cell RNAseq datasets repeatedly showing that despite the expression of the same CPM, PBMCs are highly heterogeneous. 21,22 Histone post-translational modifications with high abundance impute each other accurately We recently observed that HPTMs form conserved modules across immune cell types. 15 Therefore, we hypothesized that an HPTM could be more accurately imputed using other HPTMs than using CPMs. Indeed, each of the three ML models demonstrated higher accuracy for imputing an HPTM using other HPTMs (task 3) instead of CPMs (task 1) ( Figures 1E and S2). The mean R 2 for LR, kNN, and NP were 0.49 (range: 0.02-0.87), 0.54 (range: 0.03-0.88), and 0.58 (range: 0.09-0.90), respectively. NP had significantly higher R 2 compared to LR and kNN (p < 4e-05). Interestingly, HPTMs imputed with high R 2 were those with higher abundance and low variability across cells ( Figure 1F). It is important to note that the NP models are unaware of the mean abundance or variability of these HPTMs because the abundances of each HPTM are independently scaled by the subject before being used in the models. Taken together, our results demonstrate HPTMs with high abundances are highly predictive of other HPTMs, and strongly associated with each other.
Interpretation of neural processes models reveals histone post-translational modification association networks The higher accuracy of kNN and NP in imputing an HPTM from other HPTMs (task 3) compared to LR suggested non-linear associations between HPTMs, which are modeled better with kNN and NP than LR. To infer pairwise HPTM associations and their directionality, we developed an interpretation algorithm using noise perturbations (STAR Methods). Briefly, without re-training or modifying a kNN or NP model, we imputed an HPTM Y by systematically replacing the abundance of each input HPTM X i , one at a time, with noise. We quantified the strength of the directed association from X i to Y as the percent decrease in R 2 , with 100% denoting the strongest association and 0% denoting the weakest association. The pairs of HPTMs with strong associations formed directed networks for each panel.
Interpretation of NP models identified several known directed interactions between HPTMs in both panels.
In the acetylation panel, the strongest directional interaction was between histone H3 cleaved at threonine 22 (cleaved H3T22) and H3.3S31ph (Figures 2A and S3 iScience Article R 2 reduced by 55%. The directionality of the identified association suggests biologically relevant causality, which is in line with our recent results showing that the abundance of H3.3S31ph is impacted by cleaved H3T22 in a cell line genetically modified to prevent cleavage of H3T22. 14 We also found that perturbing H3K9ac abundance reduced the accuracy of imputing PADI4 abundance, suggesting the functional interaction between the two, which is in line with results from Kolodziej et al. demonstrating a locus-specific interaction between H3K9ac and PADI4. 23 For the methylation panel, perturbation analysis identified three distinct network components ( Figures 2B  and S3). Asymmetric and symmetric arginine dimethylation (Rme2asy/sym) strongly influenced each other and were disconnected from the lysine methylation network. This is consistent with the fact that the arginine and lysine methylations are catalyzed by different sets of enzymes. The second component consisted of H3K4me2/3, which is known to promote gene expression and occur near promoters and enhancers. The third component, on the other hand, contained known markers of gene repression and heterochromatin, including H3K27me3, H3K9me2, H4K20me3, and macroH2A. In this network component, we found that perturbing H3K9me2 abundance reduced the accuracy of imputing H3K27me3 abundance, which is in line with H3K9me2 and H3K27me3 cooperating to maintain heterochromatin. 24 In addition to these known interactions, we identified directional trans histone interactions between cleaved H3T22 and H4K16ac, and between H3K9ac, H4K5ac and H2BK5ac. We also found directional cis interactions between modifications of the same amino acid such as H3K36me1/2/3 and H4K20me2/3 iScience Article (discussion). For several HPTMs imputed with R 2 > 0.5 (e.g., H3K36me2/3, Rme2asy/sym, cleaved H3T22, H3.3S31ph, H4K16ac), there was at least one input HPTM leading to >50% decrease in R 2 when its abundance was replaced with the noise, suggesting highly conserved interactions between these HPTMs ( Figure S3). Conversely, several HPTMs imputed with R 2 > 0.5 (e.g., H3K36me1, H4K20me2/3, H3K9ac, H4K5ac, PADI4, macroH2A) had no such inputs. Since the latter HPTMs were accurately imputed, detected associations suggest interactions between their inputs (discussion). Taken together, our analyses identified several known associations between pairs of HPTMs together with their directionality. We also identified several associations that lead to new hypotheses to further explore the crosstalk of HPTMs.
In contrast, when we applied our perturbation algorithm to the kNN models, the strength of associations was significantly lower (p = 2.2e-10; Figures 2C, S4, and S5). This is expected because kNN models are non-parametric and perturbation to any input HPTM leads to the recalibration of the kNN model, which in turn can result in lower strength of associations. In contrast, NP models do not recalibrate during the perturbation analysis. Furthermore, kNN models only identified a subset of known HPTM associations (e.g., H3K4me2-H3K4me3 and Rme2sym-Rme2asy), but failed to identify the directionality of the other known associations. For instance, kNN found an equally strong bidirectional association between cleaved H3T22 and H3.3S31ph ( Figure S5) and did not identify the association between PADI4 and H3K9ac. Overall, NP more accurately imputed HPTM abundances and better identified direction of interaction between HPTMs than kNN.
Histone post-translational modification associations are conserved across subjects Each NP model is subject-specific and learns a unique embedding for each subject from the input data. However, we averaged the inferred HPTM associations across all subjects in the test set. Thus, it is possible that our overall associations are driven by a subset of subjects. Therefore, we modified our perturbation algorithm to quantify the between-subject variability in the associations (STAR Methods). Briefly, we replaced the embedding for a subject with that of a different subject and estimated the between-subject variability as the percentage change in R 2 . Variability of 0% implies the model is independent of the subject embedding and the pair-wise associations for a given NP model are conserved, whereas variability of 100% implies the associations are distinct for each subject.
The NP models predicting an HPTM from CPMs (task 1) showed the highest between-subject variability (median = 32.1%; Figure 2D). This is in line with single-cell RNA-seq studies repeatedly showing a large amount of heterogeneity in immune cells identified using CPMs. 21,22 The NP models predicting a CPM from HPTMs (task 2) had moderate between-subject variability (median = 20.1%), which is also expected since we have observed immune cell type-specific HPTM profiles that are affected by aging. 13,14 The NP models predicting an HPTM from other HPTMs (task 3) had the lowest between-subject variability (median = 4.5%; Figure S6). Overall, our results show that the associations between CPMs and HPTMs vary between subjects, presumably affected by their different environmental exposures. On the other hand, low between-subject variability in the interactions between HPTMs suggests that most of the interactions between HPTMs are highly conserved such that changes in one HPTM lead to predictable changes in the other HPTMs.
A subset of histone post-translational modifications accurately imputed several histone posttranslational modifications Based on low between-subject variability in predicting an HPTM using other HPTMs, combined with our previous observation of highly conserved HPTM correlation modules, 15 we hypothesized that a subset of HPTMs may be more informative than others. To test this hypothesis, we ranked HPTMs using their average association strength, defined as the mean strength of all the associations containing the given HPTM as the input ( Figures 3A, 3B, and S3). Cleaved H3T22 and H4K16ac in the acetylation panel, and H3K36me1 in the methylation panel were the strongest predictors of 12 other HPTMs (mean association strength >25%), suggesting an (C) Summary of hybrid NP models that use CPMs and cleaved H3T22 (acetylation panel) or H3K36me1 (methylation panel) as inputs to impute an HPTM and comparison with models in task 1 (CPMs / HPTM). Each dot represents a unique model. Lines connect models predicting the same HPTM. One-sided, paired Wilcoxon test was used to compute the significance of the improvement in R 2 due to the hybrid models (n = 36 HPTMs per task). See also Figure S7.  Figure 1F). Hence, our results suggest that the abundance of HPTM is not confounding the average strength of HPTMs associations.
Next, we investigated whether the integration of these HPTMs with CPMs would improve the imputation accuracy of other HPTMs compared to using only CPMs as in task 1. The new hybrid NP models using 11 CPMs and either cleaved H3T22 or H3K36me1 to impute HPTMs had significantly higher R 2 (mean = 0.44, range: 0.01-0.78) than the models using only CPMs (p = 9.5e-08; Figure 3C). Importantly, the hybrid NP models had 17 HPTMs with R 2 R 0.5; more than threefold increase in accurately imputed HPTMs compared to HPTMs imputed from only CPMs (task 1). Although including H4K16ac with cleaved H3T22 and 11 CPMs (mean = 0.453, range: 0-0.78) significantly increased the R 2 (p = 0.011) compared to including only cleaved H3T22, it did not increase the number of HPTMs imputed with R 2 R 0.5 ( Figure S7). These results show that a modified panel with H3K36me1 and cleaved H3T22 measured in both panels would enable imputing multiple HPTMs across panels with higher accuracy, further increasing the power of EpiTOF. This also highlights the power of our interpretation algorithm in accurately identifying associations and ranking the predictive power of each HPTM.

Neural processes models and interpretation algorithms identify influenza vaccine-associated epigenetic changes
We demonstrated that the NP models can accurately impute HPTM abundances and infer known and previously unidentified associations for healthy subjects. Next, we investigated if the same models, without any modification or re-training, could be used to impute HPTM abundances from healthy subjects affected by external stimuli such as vaccinations. We recently showed that for healthy subjects receiving the trivalent inactivated seasonal influenza vaccine (TIV), HPTM abundances are most changed 30 days after vaccination in myeloid cells, which are associated with innate memory. 8 Therefore, we evaluated the NP models imputing an HPTM from other HPTMs (task 3) and inferred pair-wise HPTM associations in a cohort of 21 healthy subjects sampled before (Day 0) and 30 days after (Day 30) receiving TIV ( Figure 4A, STAR Methods). The imputation R 2 at both these time points were strongly correlated with those from the healthy subjects in the test set (Figure 4B, R = 0.9 and 0.87, and p < 1e-12 for Day 0 and Day 30 samples, respectively), indicating that the same models can be used to impute HPTMs without significant change in imputation performance ( Figure S8). Although the inferred association strengths at both these time points were also highly correlated with those from the healthy subjects in the test split ( Figure 4C, R = 0.94, 0.91 and p < 2.2e-16 for Day 0 and Day 30 samples, respectively), some association strengths showed large deviations from the test set ( Figure S8), suggesting that most HPTM associations are conserved but a few are significantly modified due to TIV.
To identify associations that are modified due to TIV, we compared the HPTM association networks independently inferred from the Day 0 and Day 30 samples. 12 out of 16 Day 0 associations in the acetylation panel network and 17 out of 19 Day 0 associations in the methylation panel network were conserved with identical association strengths in Day 30 samples (Figures 4D and 4E). However, in the acetylation panel, the association between cleaved H3T22 and H3K14ac, H4K16ac and PADI4, H3K18ac and H3R2cit, and H2K9ac and H4K5ac were present only on Day 0 but not on Day 30 ( Figure 4D). In the methylation panel, the association between H4K20me3 and H3K27me3, and H3K9me2 and H3K9me1 was present only on Day 0 but not on Day 30 ( Figure 4E). Across both panels, all the HPTM pairs associated on Day 30 were also associated on Day 0. Taken together, the networks revealed that the associations between 6 pairs of 12 HPTMs were significantly different on Day 30 post-vaccination.
Furthermore, across both EpiTOF panels, PBMCs from the two-time points were distributed differently in the UMAP 25 space defined by these HPTMs (Figures 4F and 4H). To further quantify these distinct patterns, we clustered 26 the PBMCs (Figures 4G and 4I), computed the median abundance of the HPTMs for each cluster in each panel, and compared the sample-wise proportion of cells in a cluster across the two-time points ( Figures 4J and 4K, STAR Methods). In the acetylation panel, the cell proportions were significantly different in clusters 3 (FDR = 7.2%), 7 (FDR = 4.1%), and 12 (FDR = 5.1%) (Figures 4J and S9). Similarly, in the methylation panel, the cell proportions were significantly different in cluster 3 (FDR = 8.4%) (Figures 4K and S9). These suggest that the changes in these HPTM associations are driven by a change in cell proportions in these clusters.
We have previously shown that the overall proportions of immune cell sub-types are not significantly different between Day 30 and Day 0 samples. 8  Thus, NP models and perturbation analysis led to the identification of HPTMs whose associations changed due to TIV, and clusters defined by these HPTMs showed significant differences in cell proportions. The observation of cluster-specific cell proportion changes, but not at the overall level, suggests epigenetic reprogramming of the immune cells in these clusters, which is in line with recent studies demonstrating that memory in the innate immune system is driven by epigenetic reprogramming. 9,28 In summary, we have shown that the NP models trained using EpiTOF profiles of healthy subjects can be used to impute HPTM abundances in post-vaccinated subjects without any change in imputation accuracy. The associations between 6 pairs of HPTMs are modified due to TIV, and the clusters of cells that possibly drive these differences in associations show differences in proportions of immune cell sub-types.

DISCUSSION
We profiled 38 HPTMs and histone variants with 11 CPMs across 2 panels in more than 55 million PBMCs from 158 healthy controls across 25 experiments using EpiTOF. Using this unique dataset, we developed a computational framework for inferring directional associations between two HPTMs. As part of the development, we evaluated three machine learning methods for imputation to conclusively show that NP models impute HPTMs and CPMs with significantly higher accuracy, measured as R 2 , compared to LR and kNN models.
Increased accuracy of NP models required larger amounts of training data, computational resources, and computational time to train than LR and kNN models. However, once trained, NP models offered several advantages over LR and kNN models. First, NP models better captured associations from the underlying data than kNN models and provided new biological insights. Second, NP models were faster than kNN models when making predictions during testing. Third, the learnings from the trained NP models can be transferred to other settings using transfer learning approaches. For instance, when a new HPTM is added to the EpiTOF panel, the data and computational time required for training can be significantly reduced by using transfer learning approaches. Importantly, the NP approach can be applied to impute other multipanel data sources such as CyTOF. Taken together, the advantages of the NP model outweigh its limitations for increased resources.
We found that the CPMs and HPTMs impute each other with low to moderate R 2 . The best imputed CPM (using HPTMs as input) was CD11c (R 2 = 0.6), which is abundant only in cells from the myeloid lineage. This suggests a distinct pattern of HPTM abundances in the myeloid and lymphoid lineages. Indeed, we recently found that myeloid and lymphoid cells follow epigenetically distinct trajectories during their iScience Article differentiation from hematopoietic progenitors. 15 However, other CPMs were imputed with R 2 < 0.5. This is not surprising, since these CPMs delineate PBMCs into broad immune cell subtypes (e.g., CD4+/CD8+ T cells, classical/non-classical monocytes) whereas chromatin states, and hence HPTM abundances, which lead to mRNA expression, show substantial variability within the same cellular population. 29 We found that the HPTMs accurately imputed other HPTMs, and the perturbation-based analysis allowed inferring directionality between a pair of HPTMs. Our analysis correctly identified several known HPTM interactions. Importantly, our analyses also identified several previously unidentified associations and HPTM networks that should be further investigated in the future studies to decode the crosstalk between HPTMs. For example, we identified directed interactions between mono-, di-, and tri-methylation of lysine 36 on histone H3. Although H3K36me2 and H3K36me3 have been well studied, not much is known about H3K36me1. 30 NSD1 is known to catalyze H3K36me1/2 and SETD2 catalyzes H3K36me3 in vivo. However, it is unclear whether the substrate for the SETD2-dependent catalysis of H3K36me3 is H3K36me1, H3K36me2, or both. 31 Our perturbation analysis strongly suggests that H3K36me1, not H3K36me2, is a stronger predictor of H3K36me3, which in turn suggests that H3K36me1 is the substrate for SETD2-dependent catalysis of H3K36me3.
We identified cleaved H3T22 and H3K36me1 as strong predictors for HPTMs in the acetylation and methylation panels respectively and showed that a hybrid model using cleaved H3T22 and H3K36me1 with the CPMs significantly improves imputation accuracy for HPTMs. We note that both cleaved H3T22 and H3K36me1 are not measured in ENCODE 32 and are understudied. 30,33 Our results indicate a pair of highly abundant yet understudied HPTMs to play a central role in HPTM associations and crosstalk.
The inferred HPTM associations also revealed several accurately imputed HPTMs (e.g., H3K36me1, H4K20me2/3, H3K9ac, H4K5ac, PADI4, and macroH2A) that had no single strong predictor. These likely suggest that several predictors interact with each other such that the removal of a single predictor doesn't significantly affect the imputation of these HPTMs. Inferring these interactions would require performing perturbation analyses on several combinations of inputs for a single output. In this work, we only focus on pair-wise associations between an input and an output.

Limitations of the study
Our study has a few limitations. First, although the networks inferred from the interpretation framework provide valuable insights about associations between HPTMs, our analysis does not identify the immune cell sub-type where these associations are identified in this proof-of-concept study. As more data become available, future studies should train NP models per cell type to potentially further increase the accuracy of imputation and derive cell type-specific HPTM associations. Second, we focus only on the pair-wise associations from an input to an output. Future studies should perturb multiple combinations of inputs to quantify the effect of their interactions in the prediction of an output.
In summary, we described a machine learning-based framework to impute abundances for a subset of HPTMs with high accuracy. We found that NP-based models, which had the highest accuracy, have several advantages compared to other imputation methods. First, NP-based models capture variability at the cellular and subject levels using subject-specific functions. Second, because the NP models are based on deep learning, current trained models can be adapted using transfer learning to impute HPTMs added to the EpiTOF panel in the future. We also demonstrated an interpretation algorithm that infers valuable insights into the underlying biological processes from the data. Our methods have the potential to be applied to other single-cell frameworks such as CyTOF and scRNA-seq to infer directed associations.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: Blinding at any stage of the study 5 experiments (35 subjects) used as test set were blinded until the computational models (LR, kNN, and NP) were optimized and locked using the training and validation samples.

Sample-size estimation and statistical method of computation
Sample size estimation was not performed because our models were at single-cell resolution, where we had >55 million cells.

Inclusion and exclusion criteria of any data or subjects
In order to account for heterogeneity in healthy subjects and demonstrate that our models are generalizable, we did not have any inclusion or exclusion criteria.
All statistical analyses were performed using R and corresponding details can be found on the figure legends. The statistical significance was computed using the non-parametric Wilcoxon test, which does not have any underlying assumptions about the data distributions.

Imputation tasks, model inputs and outputs
We evaluated models imputing a CPM or an HPTM across 4 tasks: an HPTM imputed using CPMs (task 1), a CPM imputed using HPTMs (task 2), an HPTM imputed using other HPTMs on the same panel (task 3), and an HPTM imputed using CPMs and best predicting HPTM (hybrid models). Each model also uses the abundances of H3 and H4, age, and sex of the subject as inputs.

Data scaling and splitting
We scaled each CPM and HPTM by subject to have mean 0 and standard deviation 1. Age was divided by 100 and sex was one-hot encoded. To prevent overfitting and ensure generalizability, we split the data by experiments into train, validation, and test, while keeping at least 5 experiments in each split (Table S2).

NP models: Design overview
We first describe how a trained and tested model will be used for imputing an HPTM across panel and use the example of imputing H3K36me1 (Y ) from CPMs, H3, H4, age and sex (X) (task 1) for a subject s (Figure 1C). The NP model consists of two neural networks-the encoder and the imputer. The encoder processes cells from the panel where Y is measured and known. For H3K36me1, it is the methylation panel.
We denote the cells processed by the encoder as the context cells. First, the encoder encodes each context cell: cell encoding context cell i = encoderðX i ; y i Þ; where b y j represents the imputed abundance of H3K36me1 in the j th target cell.

NP models: Training and testing
For training and testing NP models, we only used cells from the panel where Y is measured. We randomly divided the cells from this panel into two equal parts and used the first as the context cells and the second as the target cells. Due to GPU memory limitations, we split cells from each subject into 20 parts randomly and processed each as a separate subject.
We computed the training loss as the sum of two terms loss = MSE subject embedding s;context cells ; subject embedding s;target cells where MSEðÞ is the mean squared error , subject embedding s;context cells˛R 512 and subject embedding s;target cells˛R 512 are computed by the encoder using the context and target cells respectively, and Y target cells˛R nt and b Y target cells˛R nt are the measured and the imputer imputed values respectively for the output corresponding to the nt target cells. The first term corresponds to the encoder network loss and the second the imputer network loss.
We used the coefficient of determination on the test split to assess the model performance: NP models: Architecture The encoder and imputer networks are fully connected with 2 hidden layers of dimension 256 (Table S3). The input dimensions for each network varies by task. The encoder output has dimension 512 and the imputer output is a scalar.
kNN and LR models kNN models impute Y for a target cell based on its nearest neighbors in the context cells. To be comparable with the NP models, the division of cells into context and target groups are the same for both NP and kNN. We evaluated k˛f1; 5; 10; 20; 50; 100; 200g on the validation split and used the best k for each model on the test split.
LR models impute Y for a target cell based on the learned linear combination of the inputs X. The LR models do not use the context cells.

Perturbation-based algorithm to infer HPTM associations
To infer the association strength between an input variable X i and output Y , we designed a perturbationbased interpretation algorithm. First, we replaced the measured values of X i in the target cells for subject s with random values drawn from N ð0; 1Þ and imputed the output variable: b Y s;X i perturbed = imputer À X s;X i perturbed ; subject embedding s Á :

(Equation 6)
We then evaluated R 2 after input perturbation: and then obtained the strength of the directed association from X i to Y : association strength X i to Y = R 2 Y ;unperturbed data À R 2 iScience Article Constructing directed HPTM association network from inferred pair-wise association strengths To obtain the association network for HPTMs in a panel from the pair-wise association strengths, we first removed HPTMs imputed with R 2 % 0:5. Among the associations between the remaining HPTMs, we calculated the threshold strength as threshold strength = meanðasscn: strengthsÞ + std:dev:ðasscn: strengthsÞ: (Equation 9) We then removed all associations with strength lesser than the threshold strength and combined the remaining associations into a single directed network.

Perturbation-based algorithm to infer subject-wise variation in NP models
To infer the subject-wise variation in an NP model, we replaced the computed subject embedding s for subject s with those from a different random subject b s and imputed the output variable: b Y s;subject perturbed = imputerðX s ; subject embedding b s Þ:

(Equation 10)
We then evaluated R 2 after input perturbation: subject wise variation Y = 0% implies the model, and the underlying associations, are independent of the subject. The greater this value, the more the variation in the underlying associations across subjects.

Dimensionality reduction and clustering for studying the effects of TIV
We computed UMAP 25 projections on a subset of 1,000,000 cells (250,000 randomly sampled cells from each panel and time point) using n neighbors = 15 and min dist = 0:1. We performed phenograph clustering 26 on the UMAP space using a subset of 100,000 cells (25,000 randomly sampled cells from those used for UMAP from each panel and time point) using k = 1000. We then transformed every cell from both the time points to the UMAP space without modifying the UMAP projection function and assigned clusters based on the 5 nearest neighbors in the subset where the clusters were computed. We computed the cluster-wise median HPTM abundances and sample-wise cell proportions using all the cells.